function [co] = replicate(co, cell_replicate)
%REPLICATE Summary of this function goes here
%   Detailed explanation goes here

co.LatticeVectors = diag(cell_replicate)*co.Lattice.LatticeVectors;


a = {1:cell_replicate(1), 1:cell_replicate(2), 1:cell_replicate(3)};
cell_replicate_indices = allcomb(a{:});
num_cells = size(cell_replicate_indices, 1);

positions = zeros(length(co.Basis.BasisVectors)*num_cells, 3);
%positions(1:co.Basis.NumAtoms, :) = co.Position;

atoms = cell(length(co.Basis.BasisVectors)*num_cells,1);

for i = 1:num_cells
  positions((i-1)*co.Basis.NumAtomsTotal+1:i*co.Basis.NumAtomsTotal,:) = ...
    co.Basis.BasisVectors*co.Lattice.LatticeVectors +...
    repmat(sum(diag(cell_replicate_indices(i,:)-1)*co.Lattice.LatticeVectors, 1), co.Basis.NumAtomsTotal, 1);
  %co.Positions*diag(cell_replicate_indices(i,:));
  atoms((i-1)*co.Basis.NumAtomsTotal+1:i*co.Basis.NumAtomsTotal) = co.Basis.Atoms;
end


[co.Positions, co.Atoms] = sort_atoms_by_type(positions, atoms);

end

